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Abstract 

We obtain the distance of closest approach of the surfaces of two arbitrary elUpsoids valid at 
any orientation and separation, measured along their inter-center vector. This directional distance 
is derived from the Elliptic Contact Function. The geometric meaning behind this approach is 
clarified. An elliptic pair potential for modelling arbitrary mixtures of elliptic particles, whether 
hard or soft, is proposed based on this distance. Comparisons to Gay-Berne potentials are discussed. 
Analytic expressions for the forces and torques acting on the elliptic particles are given. 
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I. INTRODUCTION 



Self-assembling systems of current interest in such diverse molecular electronic 

materials to biological systems often involve molecular units or supra- molecular structures 
that are highly anisotropic in shaped. Typical examples range from lipids in biological mem- 
branes^, to alkanethiols in self-assembled monolayers^., to carbon nanotubes and inorganic 
nanorodsA. In the interplay between accuracy, and, simplicity and computational efficiency, 
such fairly rigid units are often approximated as ellipsoids. 

The importance of the geometrical anisotropy in pair-wise interactions has been recog- 
nised early on. Berne and Pechukas introduced the Gaussian Overlap Potential (GOP)^ 
whose generalisation led to the widely used Gay-Berne potential (GB)^. The main idea 
behind this approach is the representation of a pair of particles by their joint elliptically 
stretched Gaussian distribution centered around the molecular centroids. The GB potential 
has been extensively used for the modelling of the phase behaviour of liquid crystals^. There 
is renewed interest in the possibility to study dynamical behaviour, as for example in lipids 
in biological membranes^, and in side-chains in coarse-grained protein force-fields^. Part of 
the appeal of the form of the GB potential is that it allows for analytic derivation of forces 
and torquesiS acting on the particles which simplifies the modelling. The drawbacks are its 
lack of generality: for example, in its original form, it is only applicable to identical uniaxial 
elliptic particles, although it has been recently extended for the special case of particles with 
biaxial symmetry^. Another extension of the GB potential was made for an "ellipsoid in 
the sea of spheres" scenarioi^, where the main semi-axis of an elliptic particle is much bigger 
than the radius of the spherical particles interacting with it. However, as the geometry of the 
particles vary, so do the number and value of the parameters introduced in the potentials. 

A different approach to this problem relies on the Elliptic Contact Function (ECF) in- 
troduced by Perram and coworker a ---'— which emphasises its geometrical aspects. The ECF 
approach calculates the distance of closest approach of two ellipsoids with given orientations. 
It gives the correct position of contact when the two ellipsoids are in tangent contact. There- 
fore it is often referred to as the "hard ellipsoids" approach. It exploits the algebra behind 
the representation of ellipsoids as quadratic forms to write the problem as an optimisation 
task that can be solved efficiently. The GOP potential is known to be closely related to the 
Elliptic Contact Potential (ECP) derived from the ECFii*i^. 
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Although apphcable to any mixture, the ECF has not yet been widely used in applications. 
The reasons may be due to the lack of simplicity of the ECF approach and the lack of clarity 
of its geometric interpretation. Furthermore, the potential derived from this approach has 
several drawbacks: it does not distinguish in energy between different relative orientations, it 
does not become isotropic at large separations and it artificially keeps the elliptic shape of the 
potential along the longer semi-axes of the ellipsoids. In this work, we clarify the geometric 
meaning of the ECF approach and the parameters associated with it. We show how the 
true distance of closest approach of the surfaces of two ellipsoids can be approximated well 
directly from the ECF by the directional distance of closest approach along their inter- 
center direction R. This allows us to develop a new type of elliptic potential applicable 
for mixtures of ellipsoids and/or spheres. We show that this can be done for any size, any 
orientation, and "hard" or "soft" particles. In all cases, the potential behaves isotropically 
at infinite separations and addresses the drawbacks discussed above. A comparison with the 
GB potentials is also presented. Finally, we derive analytical expressions for the appropriate 
forces and torques that make Molecular Dynamics (MD) simulations possible. 

In Section |n] we give the preliminaries, namely review the main aspects of the ECF 
from the literature, in Section ITTTl we show how to obtain the directional distance of closest 
approach of two ellipsoids by the value of the ECF and compare it to the GB potentials. In 
Section HVl we show how it leads to a LJ-type of elliptic potential and in Section|3we derive 
analytic expressions for the forces and torques acting on the particles due to the suggested 
potential. Finally, in Section IVII we discuss the implications of this approach and comment 
on its advantages and drawbacks. 

II. PRELIMINARIES: THE ELLIPTIC CONTACT FUNCTION (ECF) 

The shapes of two anisotropic particles "A" and "B" centred at points r and s are 
represented by the ellipsoids A{x) = 1 and B{x) = 1, respectively (Fig. H]). A{x) and B{x) 
are quadratic forms given by 

A{x) = {x — r) ■ A ■ {x — r) 

and 

B{x) = {x - s) ■ B ■ {x - s). 
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The matrices A and B can be expressed as 

i=l,2,3 i=l,2,3 

where Ui are the unit orthogonal vectors along the three semi- axes of ellipsoid "A" with 
lengths tti and, Vi are the unit orthogonal vectors along the semi-axes of ellipsoid "B" with 
lengths bi. The symbol "®" in is used to define the matrices A and B in dyadic form. 

The ECF is a measure of proximity of two ellipsoids, originally presented and sub- 
sequently reformulated iiiM. The equivalence of the two definitions is shown in Appendix 
A. We begin with the latter formulation because it allows us to obtain a clear geometrical 
interpretation. To this end, we define S{x,X) as an affine combination of quadratic forms 
A{x) and B{x), 

S{x,X) = XA{x) + {l-X)B{x), (2) 

where A is a parameter from the interval [0, 1]. The ECF is defined as a solution to the 
following optimisation problem: 

F{A, B) = max min S{x, X). (3) 

The minimum x{X) of S{x, X) for each value of the parameter A can be found from 

VS{x, A) = 2 {AA ■ (x - r) + (1 - X)B ■ {x ~ s)} = (4) 

as: 

a;(A) = {AA + (1 - X)B}-^ ■ {XA ■ r + {1 - X)B ■ s} . (5) 

As a result, the optimisation problem given by Eq|21is simplified to an unconstrained max- 
imisation of the scalar function S{x{X), A) by the scalar parameter A: 

F{A, B) = max S{x{X), A). (6) 

A 

Note that the curve x{X) : A G [0, 1] with x{0) = s and = r connects the centres r 
and s of the two ellipsoids, as can be seen in FigQ Along this curve, the gradient vectors 
VA{x) and VB{x) are parallel. We will find later the following notation useful 

XVA{x{X)) = 2XA- {x{X) -r) = X{X) , (7) 
(1 - X)VB{x{X)) = 2(1 - X)B ■ {x{X) - s) = -X(A). 
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Earlier, it was shown (— ) that the function S{x{X),X) has a unique maximum on the 
interval A G (0, 1). The value Ac at which the function reaches its maximum is called the 
contact parameter and the point = x{Xc) is called the contact point. The derivative of 
the fimction S{x{X), A) with respect to A is 

S'ix{X), A) = {Aix{X)) - Bix{X))} + 
2x'{Xf ■ {XA ■ {x{X) - r) + (1 - X)B ■ {x{X) - s)} . (8) 

The term in the second curly braces vanishes along the curve x{X) due to Eq. ((D). Hence, 
at the extremum, 

S'{x„Xc) = A{x,)-Bix,) = 0, (9) 

from which we can see that the contact point Xc is the intersection of the curve x{X) with 
the surface A{x) = B{x) (Fig. Q). EqU leads to 

Aix,) = B{x,). (10) 

Substitution of the last equation back into Eq. (j21) gives the following simple interpretation 
of the ECF value 

F{A,B) = A{x,) = B{x,). (11) 

The contact point Xc lies on the ellipsoid A{x) = F{A,B) due to Eq. (Illj). Similarly 
for B{x) the contact point x^ lies on the ellipsoid B{x) = F{A, B). The two ellipsoids are 
in tangent contact (Fig. due to Eq. (@)). As a result, the value of the ECF serves as a 
criterion of approach of any two elliptic particles. If the value F{A, B) is below unity, then 
the original ellipsoids A{x) = 1 and B{x) = 1 overlap and vice versa. When F{A, B) is 
equal to one, the original ellipsoids "A" and "B" are in tangent contact. 

The contact parameter Ac can be found numerically as a solution of Eq. (fTUI) by ap- 
proaching the surface A{x) = B{x) along the curve x{X): 

X,:A{x{X,))-B{x{X,))=0. (12) 

The convergence of this iterative process is usually fast. Numerical aspects of this iterative 
process as well as ways to speed it up will be discussed elsewhere. 

The value of the ECF can be further expressed in the form F{A,B) = R'^f{A,B), 
where R is the inter-center distance between particles "A" and "B" (See Appendix A) and 
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/ depends only on orientation, which leads to the definition of the Perram-Wertheim (PW) 
range parameter apw{A, B) (— ) as: 

(Tpw{A, B) = ^ = ^ (13) 
lF{A,B) Jf{A,B) 



There is a well-known close relationship ^^i^^ between this parameter and the Berne- Pechukas 
(BP) range parameter, aBp{A, B), first introduced in^ for the GOP. cppiA, B) can be 
expressed in terms of S{x{X), A) (see Appendix B) by substituting Ac = 1/2, 

aBp{A, B) = ^ (14) 

IS{x{l/2)M2) 



The BP range parameter is sometimes considered as a mean value approximation of the PW 
range paramete r^^i^^ , where Ac is approximated by 1/2 (see Eg. [T^. Additionally, Eq. 
is often treated as the distance of closest approach of two ellipsoids, therefore Eq. ()14|) is 
believed to be an approximation of the distance of closest approach fEq. I15j) . We will return 
to this point in Sec. 1111 CI 



III. APPROXIMATIONS AND GEOMETRICAL INTERPRETATIONS OF THE 
DISTANCE OF CLOSEST APPROACH OF TWO ELLIPSOIDS 

We now show how to obtain the distance of closest approach of the surfaces of two 
ellipsoids using the ECF value. The true distance of closest approach, (i, is a solution of the 
following minimisation problem 

d[A,B) = min\\xa — xi,\\, (15) 

subject to two constraints 

A{xa) = 1, B{A,) = 1, (16) 

namely, that Xa and x^ are points on the surface of ellipsoids "A" and "B" , respectively. 

We will now show how to approach this distance using the ECF value, both from above 
using the geometrical properties of the ECF, and from below, using a mechanical analogy. 
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A. The directional contact distance, dn 



An intersection of the line segment between the contact point Xc and the center r with 
the surface of the eUipsoid A{x) = 1 (Fig. ^ can be found using the value F{A, B) as: 

Xa = r+{x,-r)F{A,BY^'\ (17) 

The point Xa will be called a sub-contact point of A. The sub-contact point of B, Xb can be 
found in the same way, namely as the intersection of the line segment between points Xc 
and s with the surface of the ellipsoid B{x) = 1 

Xb = s+ {xc- s)F{A,B)-^/^. (18) 

The vector {xi, — Xa) between the sub-contact points of A and B is parallel to the inter-center 
vector R and is given by 

{xb - Xa) = {s-r){l- F{A, B)-i/2) = RdniA, B) , (19) 

where 

dniA, B) = R{1 - F(A, B)-^'^) =R- apw{A, B). (20) 

The geometrical meaning of the PW range parameter becomes now clear from Eq. ()20|) : it 
is the sum of projections of vectors {Xa — r) and (s — Xb) on the inter-particle vector R or 
simply the length of the vector R — [x^ — Xa) as soon as {x^ — Xa) is parallel to R (Fig. 
(121)). The distance d^ (Eq. I^UI) has a meaning of the shortest directional distance between 
two ellipsoids measured along the direction of their inter-particle vector R. 

To see this note that on the surface ^(a:;) = B{x) any affine combination S{x,X) with 
any value of the parameter A is equal to both A{x) and B{x). At the same time for A = Ac 
we have shown, that S{x, Ac) reaches it's minimum at the contact point x^.. As a result, the 
optimisation problem ©-(EI) can be reformulated as a problem of finding the minimum of 
A{x) or B{x) on the surface A{x) = B{x). 

Note that Eqs. (fT7|) and ()18|) hold not only for the contact point Xc but for any point Xc 
(see Fig. ^ on the surface A{x) — B{x) = in the following way: 

x, = r+{x,-r)S{A,)-^'\ (21) 

Xb = s + {x^- s)S{x,)-^'^ , (22) 
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where S{x) means here any affine combination of A{x) and B{x). The resulting vector 
(xfo — Xa) is also parallel to the inter-particle vector R (Fig. 12)). 

ixh-Af,) = Ril-Six,)~'^''). (23) 

The minimum of the value S{xc) is F{A, B). It is easy to show, that the value (1— iS(ic)^^^^) 
reaches its minimum together with S{xc). As a result, the distance (i/j(A, B) has a meaning 
of minimum length of the vector {x^ — Xa) parallel to R. 

If we consider the value dji{A, B) as a minimum distance along any arbitrary direction 
independently from the inter-particle vector, then its minimum among all possible directions 
naturally is the distance of closest approach of two ellipsoids d{A,B). As a result the 
following inequality holds while (i(A, B) is greater than zero: 

d{A,B)<dR{A,B). (24) 



B. The distance from below, d„ 

An estimation of the distance of Eq. (fT3j) from below is given by the distance between two 
parallel planes Ua and ah-, which are tangent to ellipsoids "A" and "B" at the sub-contact 
points Xa and x^ (Fig. Hj): 

t/„(A, B) = X, - {x, - Xa) = X, - R{1- F(A, B)-^'^) , (25) 

where Xc is the unit vector in the direction of the gradient vector (Eq. |7)) 

^ ^ ^ VA{Xa) ^ VBjx,) 

" ||X,|| ||V^(a;„)|| \\VB{x,)W 

In fact, the maximum of the distance (j25j) among all possible positions of points x^ and 
xb on the surfaces of ellipsoids "A" and "B", such that the tangent planes aa and ah are 
parallel, is equal to the true contact distance d{A, B) as long as the ellipsoids "A" and "B" 
are not overlapping. 

To show this, we consider a mechanical analogy. We consider the parallel planes and 
between two fixed rigid convex bodies "A" and "B" which are separated by some distance 
and kept apart by a constant force F. The two planes will move away from each other until 
they touch the bodies "A" and "B" at points Xa and xi, respectively (see Fig. EI). The 
reaction forces Fa and Fb from the bodies "A" and "B" will act on the planes and 
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along normal vectors of the surfaces of "A" and "B" at the points Xa and X[,. The normal 
vectors at the surfaces of "A" and "B" at Xa and x^ are parallel due to the fact that Ua and 
ctfo are parallel themselves. 

The stationary/equilibrium point will be reached when the total force Fa + Fb and the 
total torque {xb — x a) x Fb acting on the system of the planes are both equal to zero. 
The stationary point will correspond to the maximum possible distance between the planes 
because of the force F that is acting to separate them. The total force is always equal to zero 
Fa + Fb = because Fa = F and Fb = —F. The total torque can be equal to zero only 
if {xa — Xh) is parallel to the normal vectors of the surfaces of "A" and "B" at Xa and Xf,, 
which is a solution of the optimisation problem for the distance of closest approach d{A, B) 
(Eq. (US))) as well. As a result the maximum distance between planes and trapped 
between two convex bodies is equal to the distance of closest approach of the surfaces of 
these two bodies. 

This leads to 

Q<dr,{A,B)<d{A,B)<dR{A,B), (26) 

while d{A^B) > 0. All three values dn{A,B), d{A, B) and dji{A,B) are equal to zero, 
when the ellipsoids are in tangent contact 

= d^{A, B) = d{A, B) = dR{A, B). (27) 

The inequality (f^Bj) does not hold when the two ellipsoids overlap. The distance d{A,B), 
the solution of the optimisation task (fT3j) -()16 p . remains zero whenever "A" and "B" overlap. 
However, the values (i^(A, B) and dn{A, B) are below zero and characterise the overlap of 
the ellipsoids. As a result, the distance d^iA^B) can be used in both "soft" and "hard" 
potentials. 

In practise, dn{A, B) is closer to the true distance d{A, B) for small separations while 
djiiA, B) is a very good approximation at larger separations and behaves isotropically for 
infinite separations. Even though the distance d{A, B) (fT3j) can be used to build a pair 
potential between two elliptic particles^^, the usage of a pair potential based on this dis- 
tance in molecular simulations is computationally expensive because it requires a solution 
of the problem (fTH|) - p(ij] at least once at each integration step for each pair of particles. 
Although the calculation of the ECF as a proximity measure still involves the solution of 
the optimisation task (0)-® for each pair of particles on each simulation step, the resulting 
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optimisation problem is simpler and computationally more efficient than the problem ()15|) - 
(fTB|) . Additionally, we will show that the usage of the range parameter of Eq. (fT^ leads to 
simple and compact expressions for forces and torques acting on the particles. 



C. On the relation between the BP and PW range parameters 

We can now return to the relation of the BP range parameter as a mean value approxi- 
mation to the PW range parameter (-^-^i^^ ). To examine this further, let's consider the value 
of S{x, A) on the curve x{X) for different relative orientations. We first consider two iden- 
tical uniaxial elliptic particles. The GOP and GB potentials were originally introduced for 
this kind of particles^. We will show that the function iS(a;(A), A) becomes symmetric on 
the interval A G [0, 1] whenever there is a point or a line or a plane of symmetry between 
the two ellipsoids. As a result the maximum F{A, B) of the function S{x{X), A) can only 
occur at Ac = 1/2. Hence, when the quadratic forms A{x) and B{x) are symmetric, the BP 
range parameter aBp{A, B) agrees with the value of the PW range parameter apw{A, B). 
However, for asymmetric configurations, this is generally not true. 

To proceed we express the shape matrix of a spheroid as 



where 21 is the length of the ellipsoid, 2d its breadth, u a unit vector of the the main axis 

3 

and E a unit tensor given hy E = J2 ^ ^j- Let's for simplicity move the origin of the 
reference frame into a centre of symmetry c (this point is not necessarily unique for different 
types of symmetry). Vectors designating directions and positions in this reference frame are 
denoted by Greek letters. A vector ^ is transformed into its symmetric image ^ by a matrix 
n, which does not change its length but changes its directions: 



For the matrix O to be a symmetric transformation, we require the following property: 



A = l-' 



(28) 



(29) 



n n = E, 



(30) 



which implies that the determinant of fi can be 1 or —1. 
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The quadratic forms A{x) and B{x) are symmetric with respect to the centre of symmetry 
c if the centre s of B{x) is a symmetric image of the centre r of A{x), so that 

r = c + r], s = c + r], fi = O, ■ -q (31) 

,and, the matrix S is a symmetric image of the matrix A, satisfying 

B = fl A fl. (32) 

Eq. (jH^ includes in addition the symmetries of particle "A". 

The possible cases of symmetry of quadratic forms A{x) and B{x) are the point symmetry 
f2 = —E, the line symmetry, fl = 21 ^ I — E, and the plane symmetry fl = E — 2n ^ n, 
where Z is a unit vector along the axis of symmetry and n is a unit vector normal to the 
plane of symmetry. In the case of line symmetry, any point of the actual axis of symmetry 
can serve as a centre of symmetry c, while in the case of plane symmetry, any point of the 
plane of symmetry can be considered as a centre of symmetry c. 

Using Eqs. ()29|l - (j32|l . it can be shown that for any point $, and its image the values of 
the symmetric quadratic forms satisfy the following equations: 

A{^) = B{i), A{i) = l3i^). (33) 

The function S of EqElnow becomes 

5(^,A) = A^(0 + (1-A)S(0= (34) 
XB{i) + (1 - X)A{i) = Sit (1 - A)) 

and its gradient is given by 

VS{l{l-X)) = n-VS{tX). (35) 

From Eq. it follows that, if a point $, belongs to the curve x{X) (Eq. (jS))), its symmet- 
ric image belongs to the curve x{X) as well. Further from Eq. ()34|) it follows that whenever a 
symmetry described by Eqs. (j2ni)-(IS2I) is present, the function S{x{X), A) becomes symmet- 
ric on the interval A G [0, 1]. As a result, the maximum F{A, B) of the function S{x{X), A) 
can only occur at Ac = 1/2. This shows that the BP and PW parameters are equivalent for 
symmetric quadratic forms, namely for symmetric configurations of equivalent particles. 
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If we now inspect the zero level of the GB pair potential depending on the relative 
orientation of two particles 

UGBiA,B)=0, (36) 

we can see that the level (j36|) does not depend on the GB strength parameters (see Appendix 
B) and it is equivalent to 

aBp{A,B) = l, (37) 

which is in turn equivalent to 

5(a;(l/2),l/2) = 1. (38) 
In the case of point symmetry the matrix B is equal to the matrix A. 

B = n- A-n= {-E) ■ A- {-E) = A, (39) 

which corresponds to "perfectly aligned" configurations including "side-to-side" and "end-to- 
end" configurations, which are often used for the adjustment of GOP and GB potentials. (In 
this work, we use the "i-to-j" notation to refer to the different configurations. For example, 
"end-to-end" is denoted by "1-to-l", "side-to-side" by "2-to-2", "side-to-end" by "1-to- 
2".) As a result, the BP range parameter predicts the value of the PW range parameter 
correctly. However, the "l-to-2" orientation of the molecules is not symmetric and here 
the approximation Ac = 1/2 fails (Fig. E}. In fact, Eq. (jHH) will be satisfied when the 
actual ellipsoids are already separated by some distance above zero due to the fact that 
in asymmetric configurations the value iS(a:;(l/2),l/2) is always an underestimation of the 
EOF value F{A,B) (Fig. EJ. This means that the size of the molecule given by ()36|) for 
the "side-to-end" configuration (and any other configuration such that Ac 7^ 1/2) will be 
increased compared to a symmetric configuration (such as "side-to-side" and "end-to-end") 
of the same molecule and the same potential. Such an increase of the volume of the particles 
might act as an artificial " ordering force" , which will force elliptic particles into symmetric 
(ordered) configurations. 

This analysis is of course far from complete. The total pair potential also depends on 
the GOP strength parameter, so that the resulting deviations of the final shape of the 
particles are rather complex. In general, we expect the described deviation to grow with 
the elhpticity of the interacting spheroids. Note, that for unequal elliptic particles, the 
approximation Ac = 1/2 fails everywhere except for some rare cases. 
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IV. AN ELLIPTIC POTENTIAL BASED ON THE DIRECTIONAL DISTANCE, 

dR 



We now use the directional distance dii{A,B) obtained in Sec illl Al to construct a new- 
elliptic potential. Here, we concentrate on the LJ 12-6 form as an example of a radial 
potential. 

The LJ potential for two different spherical particles "A" and "B" has the following 
forrni^: 

(40) 



ULj{A,B) = 4eAB 
where 



, Pab j \ Pab , 



= (cTA + CrB)/2, eAB = \JtAA^BB , (41) 

and Pab is the inter-center distance. Eq. (PT|) represents the well known Lorentz-Berelot 
mixing rules for Van-der-Waals radii a a and as of interacting particles and their depths 
parameters caa and bbb commonly used in the LJ potential. The parameter eaa has the 
meaning of the depth of the potential well for the interaction of two identical particles of 
type "a". 

The distance df{{A, B) can be consistently compared with the inter-center distance R. 
Using Eq. (pUj) . the PW range parameter o"pi4/(A, B) can be treated as the sum of radii of 
two elliptic particles at a given relative orientation. From this point of view, apw{A, B) can 
be considered as an orientation-dependent "mixing rule" for two elliptic particles similar to 
(Jab in Eq]41l Using this analogy, the following elliptic potential can be built in Lennard- 
Jones form: 



4eo[F(A,S)-6-F(A,S)-3] , (42) 

which is the Elliptic Contact Potential (ECP) developed by Perram and co-workersi^ii^. We 
can now see that this potential is a consistent generalisation of the Lennard- Jones potential 
in the case of elliptic particles when one considers the inter-particle distance R together with 
the "mixing rule" apw{A, B). The value eo has the meaning of the potential minimum. 

However, the potential ()42j) in this form has several disadvantages. First of all, it does not 
become isotropic at large separations of the particles. Secondly, the depth of the potential 
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minimum remains the same for all relative orientations of the particles "A" and "B". And 
lastly, the shape of the potential well is not realistic because it is wider along the longer 
semi-axes of the ellipsoids than the original potential. These disadvantages prevent the 
potential of Eq. (j42j) from being extensively used in applications. Here, we would like to 
build an extension to the ECP, which preserves the good features of the ECF approach but 
frees it from the disadvantages mentioned above. 

These problems come from the " Lennard- Jones like" reduced distance {R/apw{A, B)) 
used in the ECP (Eq. 021) • This distance becomes one when the ellipsoids are in contact 
and the potential goes to zero and keeps an "elliptic shape" at large separations that causes 
both the unrealistic anisotropy of the ECP at long distances and the unrealistic shape of the 
potential well. However, the directional distance dji{A, B) along the inter-particle vector 
itself becomes isotropic at large separations. Using this distance instead, a "shifted form" 
of the potential (jl^ can be built as: 



f/(A,B) =4e, 



> 12 



dR{A,B)+mJ \dR{A,B)+ao 
( i?.-i?.F(A,B)-i/2+^()) ~ (iJ-iJF(A,B)-i/2+ao 



(43) 



where (Tq has the meaning of characteristic length and is responsible for the width of the 
potential well. The shape of the potential well is now more realistic because the distance 
(iij(A, B) is a good approximation to the distance d{A, B). The potential (j43|) has the form 
of the GB potential without a strength parameter and the PW range parameter instead of 
the GOP one^. 

The depth of the wells of the potential ()43p however still equals one in any direction. One 
way to allow for variable depth of the potential minima is to use different shape matrices 
for the attractive and the repulsive part of the potential of Eq. 



f/(Ai,A2,Si,B: 



4eo 



2) - 

,R-/?Fi(Ai,Bi)-i/2+a,J ~ ViJ-iJF2(A2,B2)-i/2+<jo 



(44) 



where the "repulsive ECF" Fi(Ai,Si) and the "attractive ECF" F2{A2,B2) are calculated 
using different shape matrices 

^1 = H (^li'^i ® Bi= J2 ^li'^i ® (45) 

i=l,2,3 i=l,2,3 
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and 

A2= (^2i'^i ^2= ^'2i'"i ® ""i ■ (46) 

i=l,2,3 i=l,2,3 

A justification for tliis approacli can be given by inspecting tlie constant levels of the 
attractive and repulsive parts of a pair potential of two complex molecules, calculated as the 
sum of attractive and repulsive parts of the LJ potentials of spherical particles. It is easy to 
see that the constant levels -LJ12 and LJq do have close but different shapes. 

Although this introduces more parameters to the problem, the total number of fitting 
parameters is still reduced compared to the GB potential. For example, the empirical 
exponents v and /i present in the GB potential (see Appendix B) are excluded. Additionally, 
the small deviation of the attractive shape matrices A2 and B2 with respect to the repulsive 
ones Ai and Bi allows us in practise to use the solution of the "attractive" ECF problem- 
Ac2 and Xc2- as the initial conditions for the "repulsive" ECF problem, which reduces the 
number of iterations required for the calculation of the latter. Besides, when calculated 
separately, the repulsive ECF problem requires a much smaller cut-off value. As a result 
both ECF problems are solved only when molecules are very close. 

To illustrate the form of the potential we consider two examples. 

A. Example 1 

As a simple test case, we consider the pair potential of two identical spheroids consisting 
of 6 LJ spherical particles in a linear array^. The resulting fitted potential and the target 
LJ potential are shown in Fig. IHl Although only the values of "side-to-side" and "end-to- 
end" potential minima are fitted, the "side-to-end" configuration is also close to the target 
potential. The "2-2" or "side-to-side" potential minimum of such potential is mm/jf/22 = 
14.883 when the LJ parameters are a = 1, e = 1, and the spheres are separated by a distance 
of 2/3. and has been set to unity. Therefore a strength parameter e = 14.883 should be used 
in simulations. 

For identical uniaxial particles there is the possibility for further reduction of the num- 
ber of the fitting parameters. The choice of " z-th" " attractive semi-axes" equal to " z-th" 
" repulsive semi-axes" 

an = bu = a2i = 62* (47) 
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leaves the corresponding " i-to-i" potential minimum equal to one as expected. Then the 
parameter eo in Eq. ()43|) has the meaning of the "i-to-i" potential minimum. As a result, 
the adjustment of the pair potential (j43p for identical spheroids requires in total 5 fitting 
parameters, while in the GB potential 8 parameters are used. 

The same reduction can be used for unequal and biaxial particles as well, but then the 
fitting of the potential to a target potential become less flexible in the representation of 
ratios of different potential minima. 

B. Example 2 

We consider a more complicated example, namely the pair interaction of two non-identical 
biaxial molecules, peropyrene 6*26-^^14 and anthracene CuHio (Fig. [7j). The peropyrene-to- 
anthracene pair potential is calculated as a sum of LJ interactions between all atoms of the 
molecules after a geometry optimisation of the isolated moleculea^li. Considering only cases 
where the "i-th" semi-axis of the particle "A" is aligned with the "j-th" semi-axis of the 
particle "B", there are in total 9 "i-to-j" configurations which need to be represented. The 
objective function for the least-squares fitting was built as a weighted sum of square differ- 
ences of the potential minimum value, the position of the root U {A, B) = and the width of 
the potential well at half-depth for each of the "i-to-j" relative orientations of the molecules. 
In this example, the relative weights were chosen to represent "aligned" configurations "1-to- 
1", "2-to-2" and "3-to-3" more accurately then other "misaligned" configurations. However, 
the particular choice of the weights can be made to suite an application in mind. The target 
potential was averaged over rotations of the molecules about their inter-center vector R for 
each "i-to-j" configuration. Fig. (jHJ shows that the potential (jl^ is able to reproduce, at 
least qualitatively, the complex interaction profile of the peropyrene-anthracene pair poten- 
tial. This is encouraging since da is more accurate at longer separations rather than at short 
but the wells around the minima are still reasonably well represented. 

V. DERIVATION OF FORCES AND TORQUES 

The dynamical evolution of the system via Molecular Dynamics (MD) simulations re- 
quires the calculation of forces and torques acting on the particles due to their interactions 
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with all other particles. The suitability of a reduced representation potential for MD sim- 
ulations is mainly defined by the computing time spent on the calculation of forces and 
torques due to that potential. A numerical estimation of the derivatives of the potential is 
computationally expensive. Additionally, approximate values of the derivatives is an extra 
source for numerical errors. The possibility of expressing those derivatives analytically is a 
big advantage for a reduced representation potential. 

In this Section we derive analytic expressions for the forces and torques acting on the 
elliptic particles due to the suggested potential (jiij) that allow for efficient MD simulations to 
be performed. The resulting formulae are simple due to the simple structure of the potential 
(j44|) and the extremum properties of the ECF value (jUj). 

The following notation will be used: 



Uab = U{Ai,A2,Bi,B2) 
X, = X{\,). 



(48) 



Given 



Uab = 4eo — ^ 



where 



G^ = iR-RFr'/^ + ao)/ao, 
the force acting on particle "A" is 



dr 



2 ^-13 dGi _ d G2 



dr 



dr 



where 



and 



dG, 1 



dr (To 

dF 



dr 



R{F-'^'-l) + ^Fr^^'^ 



"^^riAiiXfi r) — Xfi 



(49) 



(50) 



(51) 



(52) 



with Aci, the contact parameter, and Xd, the contact point of the ECF task (EqslHI) " 
with matrices Aj and Bi. The derivatives of the ECF, are given in Appendix C. 
It can be shown that 

Fb = —Fa 
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and so the force conservation equation is fulfilled. Similarly, the torque acting on particle 
"A" is 

Ta = -'-^[2G^'''F,-'/\x,, -r)x X,,- 

GfF^'^\x,2 -r)x X,2] ■ (53) 

The torque conservation equation (— ) is also fulfilled 

Ta + Tb + Rx Fb = Q. 

VI. DISCUSSION 

Anisotropic particles are of interest for interrogating fundamental questions such as 
packing2^ as well as for applications in biological and nanoscale systems. Unlike spheri- 
cally symmetric particles where the distance of closest approach is immediately given by 
their inter-center distance, this is no longer true in anisotropic systems. Obtaining the dis- 
tance of closest approach as a function of orientation for any size and shape is a difficult 
task. Additionally, the potentials must capture the geometric and energetic anisotropies of 
the system consistently. 

In this work, we have considered in detail the geometry of the ECF as an approach 
measure of ellipsoid pair particles. We have shown that the directional distance of closest 
approach B), measured along the inter-center vector, can be derived directly from the 

ECF. Like the true distance of closest approach of two ellipsoids, d{A^ B), dfi{A, B) is zero 
when the ellipsoids are in contact. Unlike d{A, B, which cannot be used in its regular form 
(ll5|l - (fTH|l to characterize any overlap of the ellipsoids, dfi{A,B) is also a measure of their 
overlap. As a result, "soft" elliptic potentials can also be built. The geometrical meaning of 
the Perram-Wertheim approach parameter apw{A, B) is clarified as the shortest directional 
distance of closest approach along the inter-particle vector R. Note that the volume of elliptic 
particles is preserved within the ECF approach. 

Since this approach is valid for any two ellipsoids, an elliptic potential is suggested which 
can be used for the modelling of any mixture of elliptic and spherical particles. Examples 
of the suggested fitted potential to pair potentials of different biaxial molecules has shown 
it is able to reproduce, at least qualitatively, complex interaction profiles. Additionally, it 
correctly reduces to isotropic interactions in both shape and magnitude at long distances. 
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For questions such as how self-assembly emerges, it may be more important to sacrifice 
quantitative accuracy at short range in order to obtain a consistent long range interaction 
that does not overestimate the range of the potential and hence impose assembly artificially. 

The structure of the potential leads to fewer fitting parameters than the GB potential. 
Analytic expressions for forces and torques acting on the particles due to the suggested 
potential are derived, which make it amenable to MD simulations. However, an iterative 
solution of the Eq. (fT^ is required for each pair of molecules on each integration step 
during the MD simulation. Although this can be done efficiently, the GB potential is com- 
putationally more effective from this point of view. On the other hand, we have shown 
in Sec. 1111 G[ that the GB potential leads to deviations of the shape and of the volume of 
interacting molecules whenever the semi-axes of the ellipsoids are different or misaligned. 
The described deviations of the volume of the elliptic particles lead to an additional artificial 
"ordering" force which might be desirable in modelling of bulk equilibrium phases, such as 
liquid crystals, but may not be as reasonable when considering dynamical properties in mix- 
tures of dissimilar nanoscale particles whose properties are strongly dependent on their size 
and shape. From this point of view, there is a "trade-off" between computational efficiency 
and accuracy in the treatment of the shape and volume of the molecules. The PW approach 
parameter and the suggested elliptic pair potential are one of possible ways to resolve this. 

Another weakness of the potential of Eq. (03)), in common with the ECP potential, has 
already been discussed in-® . For instance, in the case of the interaction of two uniaxial elliptic 
particles there are infinite number of possible "side-to-side" configurations corresponding to 
rotations of particles about their inter-particle vector R. The extremities are the "parallel" 
configuration, when the main axis of the particles are aligned and the "crossover", when 
the main axes of the ellipsoids are orthogonal. For molecules consisting of a linear array of 
spherical L J particles these two configurations obviously have different values of the potential 
minima, but they are the same in the ECP and the suggested potential as well. To take 
into account these effects one should include into the potential a dependence on the local 
curvature tensors of surfaces A{x) = const and B{x) = const at the contact point Xc or at 
sub-contact points Xa and Xf,. 

A problem in connection with this path is that the curvature of the surface of an ellipsoid 
in most cases is not the best representation of the shape of a real molecule. The fitting 
of an "i-to-j" configuration of a molecule with an elliptic potential can only reproduce the 
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characteristic length in its main directions but not the total shape. A more sophisticated 
model of molecular shape is needed to take these effects into account. The ECF approach 
in its general form of Eq. (jS)) can be extended to the case of a general convex body^^. Then 
it probably should be called a "Convex Contact Function approach". The approach keeps 
the formulas for forces and torques in almost the same form. The computational efficiency 
of the Convex Contact Function extension of the ECF depends upon the description of the 
surfaces of convex bodies. Work in this direction is in progress. 
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Appendix A. The relation between the two definitions of the ECF 

The definition (jSI)-® of the ECF which is used in this article first appeared in the 
unpublished work^^ of Perram. In this Appendix we show that it is equivalent to the original 
definition given in Ref.—. 

The inter-particle vector R = {s — r) along the curve x{X) can be expressed as 



The vector X{X) itself can be found from (j54j) as a solution to the following linear equation 



Using Eq for the scaled gradient vector X{X), the value of the quadratic form S{x{X), A) 
now becomes 



R= {x{\)-r) - {x{\)-s) = 
A-iA-i ■ X{X) + (1 - Xy'B-' ■ X{X). 



(54) 




(55) 



S{x{X), X) = XA{x{X)) + (1 - X)B{x{X)) -- 
{x{X) - rf ■ X{X) - {x{X) - sf ■ X{X) -- 
X^{X) ■ {A-^A-i + (1 - X)-'B-'} ■ X{X) 



(56) 
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The quadratic form in the expression above has the same matrix as the hnear equation ()55p. 
Substituting Eq. into Eq. (jSHl), we obtain 

Six{X), A) = X^(A) ■ {A-^A-i + (1 - X)-'B-'} ■ X{X) = 

A(l - A)^^ ■ {(1 - A)A-i + XB-'}-^ ■ RR^ . (57) 

This is the original definition of the quadratic form (j57|) given in^^ ancii^. 
Appendix B: The Gay-Berne potential 

The GB potential for identical uniaxial particles (see^,— ,— ) has the following form: 

Ugb = 4eoer(A, B ) e^,{E,, E^,R) {r]'' - r]') , (58) 

V = , (59) 

{R~aBp{A,B,R) + ao) 



where the Berne-Pechukas range parameter (Jbp{A, B, R) is given by 



aBp(A,B,R] 



^[^^■{A-' + B-'y'-R 



-1/2 



(60) 



The strength parameter e2{Ei, E2, R) has the form of the square of the range parameter 
(jUn|) calculated with different shape matrices. 

1 



e2{Eu E2,R) = - [R' ■ {El + E2}-' ■ R) . (61) 
The matrices Ei and E2 are defined as 

^^1= E ^u'^'u.^Ui, E2= ^2i'^v,®v,, (62) 

i=l,2,3 1=1,2,3 

where = 1,2,3 are unit vectors along the semi-axes of ellipsoid "A" and vectors 

Vi,i = 1,2,3 are unit vectors along the semi-axes of ellipsoid "B". The parameters en 
and e2i are responsible for the potential minima of side-to-side, side-to-end and end-to-end 
configurations. 

Using the expression ()57|1 for the quadratic form S{x{X), A), the BP range parameter can 
be found by substituting A = 1/2. 
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Appendix C: Derivatives of the ECF 



The derivative of the ECF, F{A,B) (Eq. ??), with respect to position, r, for elhpsoid 
'A" is given by 



dF{A,B) dS{x„\, 



dr 



dr 



+ S'xiXc, XcJ 



dr 



dx{Xc 
dr 



(64) 



where we have used the notation introduced in Eq. ( HHj). The derivative of Eq. IMl is 
evaluated at the contact parameter A = Ac. This greatly simplifies the expression when we 
notice that the term VxS{x{X), A) vanishes for any point on the curve x{X) due to Eq. (jHj) 
and the term iS^(a;(A), A) vanishes at the extremum point due to Eq. 0. As a result, we 
are left with 

dF] _ dS{x{X),X) 
dr dr 
Similarly, for particle "B" 



-2 Xq^Ai • {Xiyl 



(65) 



dF] 
ds 



-2(1 - Xci)Bi ■ (a?ci - s) = Xci 



(66) 



The torque T^j about an axis ej and rotation angle tpj due to the potential (jl^ is 

2 



Taj 



where 



dU 



AB 



dFi 
dFo 



dtpj 

-4eo 
4eo 



_ (IUab dFj 
.^^ dF, di,j 



(67) 



2(Jo. 



R 



_R 

2^0 



(68) 
(69) 



where Gi is given in Eq. (35) (see also Eq. (37)). The derivative dFi/dipj, similarly to Eqs 
dMll-dnSl), is 

{x,i - r) . (70) 



Aci('^ci ^) 



dipj dipj 
To calculate d A/ dipj, we write A = P ■ Aq ■ P'^, where P = P{ip) is the rotation matrix 
about the axis e by the angle ip given, in dyadic form, by 



P{ip) = e (g) e + (1 - cos{ip))E + sm{ip) ex E . 



(71) 



and Aq is a diagonal matrix of the quadratic form A{x) in the "body reference frame" 
attached to the semi-axes of the ellipsoid "A" . Using the well-known expression P^ = ex P 
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for the derivative of the rotation matrix P, the derivative dA/dipj can now be expressed as 

e, X P ■ An- + P ■ An - (e, x PV = 



Bj X P - Ao- P^ - P - Ao- P^ X Sj 



(72) 



The last equahty follows from the symmetry of the matrix Aq. Using the tensor equality 
a - {b X C) = {a X b) - C , where a and b are vectors and C is a second rank tensor, Eq. 
(ffn|l becomes 



— 2Acj(iCci r)^ - (ej x Ai) - r) — 

2Acj((a;c« t)"^ X Cj) - Ai - t) = 

{{xc^ -r) X e^) - = {Xci x {x^, - r)) - Bj 



(73) 



where the last operation is just a cyclic permutation of vectors in the calculation of the 
volume of a parallelogram built by the three vectors {x^ — r), Bj and X^. 

Substituting Eq. 73 into Eq. 67, the total torque vector Ta acting on the particle "A", 

3 



- E 



i=l 



dUAB 



E {X^ X {Xci - r)) • BjBj 



- E ^ [{X,^ X {x,,. - r)) - E)] 



i=l 



E ^ [(^c. - r) X X, 



1=1 



(74) 



where E is the unit tensor E = E Cg> Bj. 
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Fig.l: The ellipsoid particles "A" and "B" , given by A{x) — 1 and B{x) — 1, are centered 
on r and s respectively. We consider the local reference frame attached to the centre of 
ellipsoid "A" ahgned along its semi-axis. We illustrate here, without loss of generality, the 
concepts behind the EOF approach in the 2D case. The original ellipsoids are scaled up 
(or down) by ^Jf(A, B) until they touch each other tangentially at the contact point Xc- 
The contact point Xc is the intersection of the curve x{\) (green line) with the surface 
A{x) — B{x) (black line). The value of the ECF is illustrated by the scaled ellipsoids 
A(xc) = F(A, B) and B{xc) = F{A, B). 

Fig. 2: The directional distance of closest approach, dR{A,B), is the minimum distance 
between the surfaces of two ellipsoids that is parallel to the inter-centcr vector R. It is given 
by the sub-contact points Xa and Xb, which are on the surface of ellipsoids "A" and "B", 
respectively. It is a good approximation from above to the true contact distance. The 
value of the PW approach parameter apw{A, B) can be understood by the geometry of the 
triangle of points r, s and Xc, as the length of the vector R — [xa- Xb) as soon as {xa- Xb) 
becomes parallel to R 

Fig. 3: The distance dn{A, B) is defined as the distance between two parallel planes 
and Ub which are tangent at any point Xa and Xb on the surface of ellipsoids "A" and 
"B", respectively. It approaches the true distance of closest approach from below. The 
maximum distance dn{A, B) coincides with the true distance of closest approach. This can 
be understood from a mechanical point of view: if the planes are kept apart by a constant 
force F, then the equilibrium point reached corresponds to the maximum distance between 
planes, which is also the distance of closest approach. 

Fig.4: A closer look at the sub-contact points Xa and tc^, which define the directional 
distance of closest approach dR. The planes and are tangent to the ellipsoids A{x) = 1 
and B{x) — 1 respectively at these points (shown with purple lines). The distance o?„ 
approaches the true distance of closest approach from below, while du approaches it from 
above. 

Fig.5: S{x{X),X) is shown as a function of the parameter A for the "2-to-2" ("side- 
to-side"), "1-to-l" ("end-to-end") and "l-to-2" ("side-to-end") relative orientations of two 
identical ellipsoids. The ECF value, Fij{A, B), is the maximum of the function S and gives 
the contact parameter Ac for each orientation. The approximation of Ac by 1/2 holds for 
symmetric configurations of identical eUipsoids such as "2-to-2" and "1-to-l" , but fails in the 
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asymmetric "l-to-2" configuration. This is expected to hold for asymmetric configurations 
in general, which includes any configuration of non-identical particles. 

Fig. 6: The resulting fitted elliptic potential of Eq. |^ (red solid curves -) and the target 
atomistic LJ potential (blue dashed curves - - ) of two identical spheroids. The original 
system consists of 6 identical spherical LJ particles separated by a distance of 2/3 and LJ 
parameters a = 1, e = 1. The "2-to-2" or "side-to-side" minimum of such potential is 
minBU22 = 14.883 which has been set to unity. Although only the values of "1-to-l" and 
"2-to-2" potential minima are fitted, the "l-to-2" or "side-to-end" configuration is also close 
to the target potential. The width parameter ctq is equal to unity. The semi-axes of the 
repulsive matrices Ai and Bi are an = bu = 0.475, ai2 = fei2 = 2.025, while the semi-axes 
of the attractive matrices A2 and B2 are given by 021 = &21 = 0.475, 022 = ^22 = 1-875. 

Fig. 7: The test molecules of the unequal biaxial case: peropyrene, C26-f^i4, and an- 
thracene, Ci/iHiQ. 

Fig. 8: "i-to-j" potential minima of the peropyrene-anthracene LJ atomic pair potential 
(red dashed lines — ) fitted with the suggested potential of Eq. EH (blue solid curves-). 
The fit was done with a weighted sum of squares difference which in this case gave more 
weight to the accuracy of the "aligned" configurations "1-to-l" ,"2-to-2" and "3-to-3". This 
choice would depend on the particular application. The following LJ parameters were used 
for hydrogen, H, and carbon, C,: cthh = 2.4A, e/^j;/ = 0.02 kcal/mol, acc = 3.4A, ecc = 
0.15 kcal/mol. The resulting repulsive shape matrix elements for Ai and Bi are an = 
7.024, ai2 = 4.006, ai3 = 1.582 and bu = 4.846, 612 = 2.841, 613 = 1.511, while the attractive 
shape matrix elements for A2 and B2 are given by 021 = 6.528, 022 = 3.829, 023 = 1.702, 621 = 
4.267,^22 = 2.405,^23 = 1-394, with width parameter ctq = 2.921 and depth parameter 
e = 25.543. 
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FIG. 1: The ellipsoid particles "A" and "B", given by A{x) = 1 and B{x) = 1, are centered on r 
and s respectively. We consider the local reference frame attached to the centre of ellipsoid "A" 
aligned along its semi-axis. We illustrate here, without loss of generality, the concepts behind the 
ECF approach in the 2D case. The original ellipsoids are scaled up (or down) by y/F{A, B) until 
they touch each other tangentially at the contact point Xc- The contact point Xc is the intersection 
of the curve x{X) (green line) with the surface A{x) = B{x) (black line). The value of the ECF is 
illustrated by the scaled ellipsoids A{xc) = F{A,B) and B{xc) = F{A,B). 
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FIG. 2: The directional distance of closest approach, dR{A,B), is the minimum distance between 
the surfaces of two ellipsoids that is parallel to the inter-center vector R. It is given by the sub- 
contact points Xa and Xb, which are on the surface of ellipsoids "A" and "B", respectively. It is 
a good approximation from above to the true contact distance. The value of the PW approach 
parameter apw{A, B) can be understood by the geometry of the triangle of points r, s and Xc, 
as the length of the vector R — {xa- x^) as soon as {xa- Xf,) becomes parallel to R 
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FIG. 3: The distance dn{A,B) is defined as the distance between two parallel planes aa and Ob 
which are tangent at any point Xa and Xb on the surface of ellipsoids "A" and "B", respectively. 
It approaches the true distance of closest approach from below. The maximum distance dn{A,B) 
coincides with the true distance of closest approach. This can be understood from a mechanical 
point of view: if the planes are kept apart by a constant force F, then the equilibrium point 
reached corresponds to the maximum distance between planes, which is also the distance of closest 
approach. 
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FIG. 4: A closer look at the sub-contact points Xa and Xf,, which define the directional distance of 
closest approach dji. The planes and Of, are tangent to the ellipsoids A{x) = 1 and B{x) = 1 
respectively at these points (shown with purple lines). The distance d„ approaches the true distance 
of closest approach from below, while approaches it from above. 
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FIG. 5: S{x{X), X) is shown as a function of the parameter A for the "2-to-2" ("side-to-side"), "1- 
to-1" ("end-to-end") and "l-to-2" ("side-to-end") relative orientations of two identical ellipsoids. 
The EOF value, Fij(A, B), is the maximum of the function 5 and gives the contact parameter Ac for 
each orientation. The approximation of Ac by 1/2 holds for symmetric configurations of identical 
ellipsoids such as "2-to-2" and "1-to-l", but fails in the asymmetric "l-to-2" configuration. This 
is expected to hold for asymmetric configurations in general, which includes any configuration of 
non-identical particles. 
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FIG. 6: The resulting fitted elliptic potential of Eq.^](red solid curves -) and the target atomistic 
LJ potential (blue dashed curves - - ) of two identical spheroids. The original system consists of 6 
identical spherical LJ particles separated by a distance of 2/3 and LJ parameters a = 1, e = 1. The 
"2-to-2" or "side-to-side" minimum of such potential is mm^?722 = 14.883 which has been set to 
unity. Although only the values of "1-to-l" and "2-to-2" potential minima are fitted, the "l-to-2" or 
"side-to-end" configuration is also close to the target potential. The width parameter do is equal to 
unity. The semi-axes of the repulsive matrices Ai and Bi are an = bn = 0.475, 012 = 612 = 2.025, 
while the semi-axes of the attractive matrices A2 and B2 are given by 021 = 621 = 0.475, 022 = 
622 = 1.875. 
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FIG. 7: The test molecules of the unequal biaxial case: peropyrene, C26-ffi4, and anthracene, 
Ci4-ff lo- 
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FIG. 8: "i-to-j" potential minima of the peropyrene-anthracene LJ atomic pair potential (red 
dashed lines — ) fitted with the suggested potential of Eq.E](blue solid curves-). The fit was done 
with a weighted sum of squares difference which in this case gave more weight to the accuracy of the 
"aligned" configurations "1-to-l" ,"2-to-2" and "3-to-3". This choice would depend on the particular 
application. The following LJ parameters were used for hydrogen, H, and carbon, C,: auH = 
2AA ■i^HH = 0.02 kcal/mol, acc = "i-^-^i^cc = 0.15 kcal/mol. The resulting repulsive shape 
matrix elements for Ai and Bi are an = 7.024, ai2 = 4.006, ai3 = 1.582 and bn = 4.846, 612 = 
2.841, 613 = 1.511, while the attractive shape matrix elements for A2 and B2 are given by 021 = 
6.528, a22 = 3.829,023 = 1.702,^21 = 4.267,^22 = 2.405,623 = 1-394, with width parameter uq = 
2.921 and depth parameter e = 25.543. 
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